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Spermatozoa flagella are known to synchronize when swimming in close proximity. We use a model 
consisting ol two-dimensional sheets propagating transverse waves of displacement to demonstrate 
that fluid forces lead to such synchronization passively. Using two distinct asymptotic descriptions 
(small amplitude and long wavelength), we derive the synchronizing dynamics analytically for ar- 
bitrarily shaped waveforms in Newtonian fluids, and show that phase locking will always occur for 
sufficiently asymmetric shapes. We characterize the effect of the geometry of the waveforms and the 
separation between the swimmers on the synchronizing dynamics, the final stable conformations, 
and the energy dissipated by the cells. For two closely-swimming cells, synchronization always oc- 
curs at the in-phase or opposite-phase conformation, depending solely on the geometry of the cells. 
In contrast, the work done by the swimmers is always minimized at the in-phase conformation. As 
the swimmers get further apart, additional fixed points arise at intermediate values of the relative 
phase. In addition, computations for large-amplitude waves using the boundary integral method 
reveal that the two asymptotic limits capture all the relevant physics of the problem. Our results 
provide a theoretical framework to address other hydrodynamic interactions phenomena relevant to 
populations of self-propelled organisms. 

I. INTRODUCTION 



An often observed yet surprising physical phenomenon is the synchronization of the pendulums of grandfather 
clocks. When two such clocks are located in close proximity, forces transmitted through a medium connecting the two 
clocks can lead to their beating in perfect synchrony [ljl. Similar synchronization can easily be obtained at home using 
two connected metronomes, with spectacular results 3|. Still more fascinating is the many examples of synchrony 
which occur in the natural world, from pacemaker cells in a heart Q , to synchronously flashing fireflies Q] . 

One particularly interesting example of synchronization occurring in nature is the observed phase locking of the 
flagella of swimming eukaryotes such as spermatozoa . These cells, typically tens of microns long, actuate slender 
flagella beating periodically in order to propel themselves in viscous fluids [9l4llj. As illustrated in Fig.[T]in the case of 
two bull spermatozoa, when two such cells swim in close proximity, their flagella are often observed to beat in synchrony 




FIG. 1: (A to C) Time-sequence showing the synchronization of two swimming bull spermatozoa. Scale bar is 25fim. Repro- 
duced from Woolley et al. 0| with permission. 
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- so much so that in Fig. [Tp the two flagella cannot even be distinguished This synchronization is biologically 
significant because it is observed to lead to an increased swimming speed for the co-moving cells, thereby providing 
a competitive advantage over cells which are not synchronized 0, • This behavior can arise purely passively, as is 
the case with the pendulums, but here the medium transmitting the forces is the fluid between the cells. While large 
systems of many bodies may be too complicated to address rigorously, and idealizations such as the Kuramoto model 
must be employed [13l |. in this paper we consider the simple case of a pair of co-swimming two-dimensional cells. We 
show that the coupling between the bodies can be obtained directly by solving analytically and numerically the field 
equations for the surrounding flow and find the occurrence of passive hydrodynamic synchronization for all but the 
most symmetric flagellar waveforms. 

G.I. Taylor first studied synchronizing flagellated cells by modeling them as two dimensional sheets propagating 
sinusoidal waves of transverse displacement [14[. With this model, he found that, for a given swimming gait, swimming 
in-phase synchronously is the conformation in which the cells swim while doing the least amount of work against the 
surrounding fluid. Left open was the question of whether the synchronization would occur passively from a random 
initial phase shift between co-swimming cells. Subsequent numerical works using an immersed boundary method and 
multiparticle collision dynamics seem to indicate that indeed synchronization could occur due to hydrodynamic forces 
alone @,0,[l6|. 



The phase locking of flagellated microorganisms is closely related to another important observed synchrony in nature, 
that of eukaryotic cilia. Cilia are short flagella typically lining the surface of a larger body and are found to beat 
in unison with a small constant phase difference giving rise to a collective motion described as metachronal waves 
[lo| . This motion provides various biological functionality including fluid transport and locomotion 17 1. Several 
models with varying complexity ha ve i ndi cated that the synchronization which manifests as metachronal beating 
can occur due to fluid forces alone 
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although, since individual cilia are not free-swimming but are attached 
to a substrate, synchronization can only occur with a load-dependent force generation. Similarly to cilia, there is 
an observed synchronization of the pairs of flagella used for propulsion on the alga Chlamydomonas [22| . Beyond 
eukaryotic flagella and cilia, hydrodynamic interactions in bacterial flagella lead to the creation of flagella bundles 
propelling the cells forward as they swim, as well as the disruption of such bundles when the cells change their 
swimming direction 
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In this paper we return to the two dimensional model first proposed by Taylor (detailed in Sec. HTJ, to describe 
the phase locking of swimming flagellated cells. The simplicity of such a model allows one to address the problem 
analytically, to extract the relevant properties that such waves must possess in order to give rise to synchronization, 
and to determine precisely what states of dynamic equilibrium will occur. We first present geometrical arguments 
which show that Taylor's purely sinusoidal sheet cannot dynamically synchronize due to an excess of symmetry 
which, when coupled with the kinematic reversibility of the Stokes equations, prevents any relative motion between 
free-swimming cells (Sec lIII[) . Real flagella possess a front-back asymmetry and we show that this feature leads to the 
occurrence of synchronization. We accomplish this by allowing the sheets to pass completely general waveforms in 
our model. We then solve the problem analytically for two asymptotic limits, first when the amplitude of the waves 
is much smaller than their wavelength (Taylor's limit, Sec. IIV|) . and then when the mean distance between the waves 
is much smaller than the wavelength (lubrication limit, Sec. W\ . We also solve the problem numerically using the 
boundary integral formulation of the Stokes equations to demonstrate the validity of the analytic formulae and to 
address the synchronization of large-amplitude waveforms (Sec. IVIj) . 

Our results show precisely how the geometry of the waveforms governs the synchronizing dynamics of the system 
(Sec. IVIIj) . We obtain simple formulae that dictate the time-evolution of the phase and the energy dissipation, and 
which indicate that while swimming in-phase results in a minimum of viscous dissipation it does not necessarily 
coincide with an equilibrium state, and indeed a dynamically stable state may maximize energy dissipation. In 
addition to the geometry of the waveforms, we demonstrate the importance of the separation of the sheets on the 
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FIG. 2: System of parallel and identical two-dimensional infinite sheets in a frame moving with the lower sheet. The sheets 
are separated vertically by a mean distance h. The top sheet, behind the bottom sheet by a phase <f> as measured along the z 
axis, moves to the right with a relative velocity Ua- 

dynamics of the system. We show that the stable conformations (and the number of them) may change with the 
distance between the cells. Notably swimming cells with front-back asymmetry are shown synchronize into either a 
stable in-phase or opposite-phase conformation when in close proximity, while some cells when further apart are shown 
to synchronize with a fixed finite phase difference, reminiscent of ciliary phase locking. A discussion and summary of 
these results is offered in Sec. IVIIII 



II. SETUP 

Our system, as illustrated in Fig. [21 consists of two parallel and identical infinite sheets, which we will call swimmers, 
separated by a mean distance h. The sheets both propagate waves of transverse displacement in the positive z direction, 
with amplitude a and speed c = w/fc, where u is the wave frequency and k is the wavenumber, and have an initial 
phase difference (j>o = kAzo (denoted positive when the bottom sheet is shifted by <fio along the positive z direction 
with respect to the top sheet). By passing these waves the swimmers propel themselves in the — z direction [ijj]. We 
consider the frame of reference moving with the bottom sheet, at speed U, and write the relative speed of the top 
sheet in the z direction as Ua- 

The instantaneous positions of the bottom (t/i) and top (2/2) sheets are thus given by 

2/1 = ag{k[z - ct]j, (1) 



2/2 = h + ag k 



z-ct + Az - / U A (t')dt' 







(2) 



where g a function describing the arbitrary waveform of the swimmers, and z is the axial coordinate in a frame moving 
with the lower sheet. We use the following dimensionless variables z* = zk, t* = tui, u* = u/c, v = v/ec, with the 
ratio of the amplitude of the waves to their wavelength given by e = ak. For convenience we use the wave variable 
x* = z* — t* and the instantaneous phase difference 4> = 4>a — k J* U A (t')dt'. Consequently the positions of the sheets 
in the moving frame are given simply by 

2/1* = eg(x*), (3) 
y* 2 = h* + eg(x* + 0), (4) 

where the arbitrary 27r-periodic function g can be written using Fourier series as 

00 00 
g(x*) = a n cos(nx*) + ^ f3 n sin(nx*). (5) 

7i — l n—1 
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Since we are concerned with the synchronization of microorganisms, we are in a low Reynolds number regime 
(Re ~ 10 -4 for the bull spermatozoa in Fig. [1]) where the fluid between the sheets is inertia free, and thus mechanical 
equilibrium for the stress tensor, er* , is written as V ■ <t* = 0. Assuming an incompressible Newtonian flow we obtain 
the Stokes equations for the dimensionless velocity field, u* = (u*,v*), and dynamic pressure, p* = pe 2 /fioj, as 

V 2 u* = Vp*, (6) 
Vu* = 0. (7) 

Physically, if the sheets are not permitted to move relative to one another, i.e. if we set f/£ = 0, then there may 
arise a horizontal hydrodynamic force f x acting on the swimmers. Conversely, if we let the sheets move freely under 
the constraint that they are force free then there may be a nonzero evolution of the phase in time, given geometrically 
as = —d(j>/dt. These two problems are of course related, as we will see, by the mobility M, as Ua = Mf x . In the 
case of a purely sinusoidal swimmer (i.e. /?i = 1, /3 n = n > 1 and a n = Vn), G.I. Taylor [l4| derived the swimming 
speed of a single sheet (the outer problem) and obtained 

U* = -r 2 ( 1 -i? e2 ) +0(e6) - (8) 

In the rest of the paper we drop the * notation for convenience. 



III. SYMMETRY 

Before calculating the hydrodynamic forces between the swimmers, it is insightful to first consider the various 
symmetry properties of the problem, and their consequences on force generation and synchronization. 

Suppose first that we have two swimmers, g\ and g 2 , whose shapes are such that g 2 is obtained from gi by a vertical 
axis reflection plus a horizontal axis reflection and a phase shift 9 (which depends on the location of the vertical 
axis), i.e. g 2 {x) = —gx(—x + 9). In that case there can be no horizontal hydrodynamic force acting between the 
swimmers, and f x = 0. To prove this result, let us assume that a force f acts on the top sheet with Ua = (since 
V • a = the force on the bottom sheet must be equal and opposite in sign) . We then perform a reflection of the 
entire conformation about the vertical axis then horizontal axis, followed by a reversal of the kinematics (see Fig. [3] 
for an example). The resulting system is identical except the sign of the force has reversed, f — > — f, a contradiction 
unless f = (then f x = 0). In particular, if the sheets are identical, then there can be no synchronization if the 
identical shapes of the waveforms satisfy g(x) = —g(—x + 6). A subset of these shapes are sheets that are invariant 
under both vertical axis reflection g(x) = g(—x + 9) and horizontal axis reflection g(x) = —g(x + 7r); the simplest 
example of such shape is a pure sinewave (fi\ = 1, /3„ = n > 1 and a n — Vn), which is Taylor's original geometry 
[14I fig . Since such an arrangement has both vertical and horizontal axis symmetry it will not passively synchronize 
in a Stokesian flow [2^ . Similar excessive geometrical symmetries have also been observed to curb any phase- locking 
in other swimmer models 27 



A further generalization of the argument may be obtained by noting that in two dimensions the outer problem 
can balance no force and hence each side of the swimmer must be force free. This decoupling of the inner and outer 
problem means that it is only the fluid between the two sheets that drives synchronization, if any. Thus if two 
swimmers do not phase-lock, a similar arrangement of more than two swimmers will not either - a result that cannot 
be obtained by symmetry alone. 

In order to possibly obtain a passive synchronization between the swimmers we must therefore either (1) have a 
geometry such that g{x) ^ —g(—x + 9), or (2) remove the kinematic reversibility of the flow equations. Since we 
are considering here microorganisms in a Newtonian fluid, the latter is a property of the problem that we cannot 
escape. If our model were at finite Reynolds number, or in a viscoclastic fluid, then this constraint would naturally 
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FIG. 3: A system of two identical and parallel swimmers which has a stabilizing force (top left) becomes destabilizing (bottom 
left), under two reflections - first about the vertical axis (R v ) then about the horizontal axis (Rh) - combined with an application 
of kinematic reversibility (KR), yet the boundary conditions remain identical, hence the force must be zero. Symmetric 
waveforms can thus not synchronize. 



be removed and symmetric swimmers could synchronize [30(. In a Stokesian flow we must thus have a geometrical 
asymmetry. 

Most swimming microorganisms, such as spermatozoa, possess a cell body and thus have a very natural front-back 
asymmetry. In addition, some spermatozoa pass waves along their flagella which increase in amplitude from head to 



tail, leading to another type of front-back asymmetry |31|. In contrast, swimmers whose flagellar waveforms or body 
is asymmetric with respect to the horizontal axis experience viscous torques, and thus cannot swim straight. It is 
therefore natural for us to focus on waveforms which are symmetric about the horizontal axis, but not the vertical. 
As a result of this horizontal axis symmetry, the horizontal component of a force between the swimmers must be an 
odd function of the phase <j>, fx{—<t>) = — fx(<l>)> and thus there must always be a fixed point at <f> = 0, i.e. /(0) = 0. 
In addition, because the force is 27r-periodic, then cf> = tt must also be another fixed point, i.e. f(jr) = 0. 

As a side note, we observe that because of the kinematic reversibility of the Stokes equations, a change in the 
direction of wave propagation yields a reversal of forces f — > — f. Reversing the direction of wave propagation is 
geometrically equivalent to reversing the front-back asymmetry of the waveforms which must therefore also reverse 
the forces on the swimmer. 

In order to gain physical intuition in the synchronization process, we now characterize the force generation and 
subsequent synchronization between the two sheets analytically by focusing on two asymptotic limits. We first consider 
in Sec. IIVI the limit in which the amplitude of the traveling waves is much smaller than their wavelength. The limit 
in which the distance between the swimmers is much smaller than their wavelength will then be considered in Sec. [V] 
Additionally we solve the problem numerically using the boundary integral formulation of the Stokes equations in 
Sec. I VII to validate our asymptotics and address large- amplitude swimming. 



IV. SMALL AMPLITUDE EXPANSION 

Because the model is two dimensional we may introduce the stream function formulation and write u = 
{dip/dy, —dip/dx}. In this manner the continuity equation is automatically satisfied and the Stokes equations reduce 
to a biharmonic equation in the stream function 



vV = o. 



(9) 
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We assume in this section that the amplitude of the traveling wave is much smaller than their wavelength, t « 1, 
and look to solve this problem by seeking a regular perturbation expansion in powers of e, tp = ^e™V , n- Because 
of the symmetry of the problem there is no difference in the boundary conditions if we change e — > — e as this is 
equivalent to taking x -¥ x + it, this then naturally precludes the possibility of a synchronizing force appearing at all 
odd powers in e. 



A. Boundary Conditions 

We wish to prescribe a wave of transverse displacement to each sheet. However, doing so requires the material 
composing the sheets to be extensible as material points will accelerate relative one another. If we wish to forbid this 
relative motion, we may require the sheet to pass waves in an inextensible fashion. This may be visualized as material 
points moving along a conveyor of static shape, when observing the sheet in the wave frame 32 



For extensible sheets the boundary conditions are given simply by the time derivatives of the waveforms namely 

u \ y = Vl = 0, (10a) 
v \ y=Vl = dyi/dt, (10b) 

U \y = y 2 = 0, (IOC) 

v \ y = V2 = dy 2 /dt. (lOd) 

For inextensible inextensible sheets the boundary conditions are given by 

u\ y=yi = 1 Sr. cos 6 \y=y x ) (11a) 

v\ y=yi = -S R sin0 \ y=yi , (lib) 

u \y=V2 = 1 _ Srcos0 \ v =y 2 +Ua, (11c) 

v \ y= y 2 = -S R sin0 \y = y 2 . (lid) 

where the angle, 9, is defined by tan# = dy/dx hence 

cos 9 = ■ 1 (12) 
sin# = y cos 9, (13) 

and the material velocity (in the wave frame), Sr., is ratio of the length of the sheet to its wavelength multiplied by 
the wave speed, or 

1 f27T / / a x 2 

(14) 



B. Expansion 

Since we know that an expansion can yield a synchronizing force only at even powers in amplitude one would hope 
to see a relative force generated at order e 2 . We actually show below that for any waveform g{x), the force is zero 
at order e 2 , and hence a perturbation expansion must carried out to order e 4 in order to obtain the synchronizing 
dynamics. 
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The governing equation at 0(e) is 



and the boundary conditions are given by 



1. Flow at O(e) 



v-0i | y =o = Vg(x), 

W-i \ y= h = S7g(x + <j))+e y U A1 , 



(15) 



(16a) 
(16b) 



where e y denotes the unit vector in the y direction. We note that the boundary conditions at 0(e) are the same for 
both extensible and inextensible motion. 

The biharmonic equation can be solved by repeated separation of variables. The general solution may be expressed 



ipx(x, y) =Ai, + B lt0 y + C lfl y 2 + £>i, y 3 + (Ei, + F 1>0 y + d.oy 2 + H lfi y 3 )x 

OO 

+ ^2 [( Al > n + B i^y) sinh(ny) + (Ci >n + Di >n y) cosh(ny) cos(nx) 
n=l 

OO 

+ ^ [(-Ei.n + Fi,nV) sinh(ny) + (Gi,„ + H^ n y) cosh(ny) sin(raa;), 



(17) 



where for the constants A through H, the first subscript refers to the order in the expansion (here, first order) and 
the second refers to the corresponding Fourier mode. We can immediately discard the terms linear in x due to the 
periodicity of the problem. 

From the first order boundary conditions, Eq. (|16[) . we get that the solution to the biharmonic equation may be 
written analytically as 



where 



ipi = ai } o(y) + ^2 \a>i,n(v) cos(nx) + bi, n (y) sin(na;) 

n=l 

(U A1 - m^oh 2 ) y 



and 



ai,o(v) 

h, n (y) 

Pn{y) = 



2h 



2P n (y)(a n cos(n4>) +/3„sin(n0)) + a n Q n (y), 
2P n (y){j3 n cos(n0) - a n sin(n0)) + PnQniv), 



n 2 hy cosh(nh) + smb(nh)ny 



(18) 



(19) 
(20) 
(21) 



cosh(nw) 

2n 2 h 2 -2smh 2 (nh) J V ' 
(l + hn 2 y) smh(nh) + hncosb(nh) 



.GO 



2n 2 h 2 -2 smb. 2 (nh) 
2nh + 2ny sinh 2 (nh) + smb(2nh) 



sinh(ny), 



(22) 



1 



2n 2 h 2 -2 smb. 2 (nh) 
2hn 2 y + ny smb(2nh) 



sinh(ny) 



cosb(ny). 



(23) 



2n 2 h 2 -2 sinh 2 (nh) 

In order to solve for the unknown constant Di q we turn to dynamical considerations. To compute the force on the 
bottom sheet we note that we are free to move the integral along the surface of the sheet S, to any surface parallel 
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to the x axis. This can be shown by integrating V • cr — over the area between the sheet and any such surface and 
using the periodicity of the problem. Alternatively, this can be shown by expanding as follows 



fx \y=yi 



/ (cr • n) \ y=yi dS 
Js 

(^xy £Q (jK)o~xx\ \y—yi 



dx 



_n— 1 „/™\n 



Qyn-l 



,=o dx 



2tt 



a 



xy ly=0 



=0 dx. 



We will use the result given by Eq. (|24|) repeatedly throughout the paper 
Using the above we find that the force on the bottom sheet to O(e) is 

pihr / d 2 ^ 1 d 2 tp 



fix \y=yi 



r ) \y=o dx 



o y dy 2 dx 2 
= 2tt< )0 (0) 

2tt 

= -6-Dio/wr + -Y- Uai, 
n 

while the force on the upper sheet is similarly 

fix \ y =y 2 = — 27ra'/ (/i) 

= -6D 10 lnr ~ -j-Uai- 
h 

Hence we see that only the second derivative of the zeroth Fourier mode contributes to the force. 
Finally, integrating mechanical equilibrium, V • cr = 0, between the two sheets leads to the equality 



(24a) 
(24b) 
(24c) 
(24d) 



(25a) 
(25b) 
(25c) 

(26a) 
(26b) 



(27) 



f \y=yi +f \y=y2 

where f = J" cr ■ ndS. Taking the ^-component we find f\ x \ y=yi = - 
we must have D\ q = 0. The force on the top sheet is therefore 

2tt 

fix = —-^Uai- 
h 

If Uai = then there is no phase locking force. Conversely if the sheets are force free then Uai = 0. There is thus 
no synchronization at 0(e), as expected from the e — ► — e symmetry. 



-fix \y=yi an d in order to satisfy this relationship 



(28) 



2. Flow at 0(e 2 ) 



The governing equation at 0(e 2 ) is 



while the boundary conditions are given by 



VVa = 0, 



d^i 



(29) 



vv> 2 | w =o = -g(x)V [ ) 



2 \y=h 



I \'(x) 2 dx 
27T Jo 

e y U A 2-g(x + <P)w(^j \ y=h 



2tt 



2tt 



'(x + cf>) 2 - — / g'(x + <j>) 2 dx 



(30a) 



(30b) 
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where the terms in the square brackets represent the contribution from the inextcnsibility constraint. 
We may write the solution to the biharmonic equation as 



V^2 = a 2 fi{y) + ^ fl 2,«(y) cos(nx) + ^ b 2 , n (y) sin(nx). 



n=l 

The zeroth mode is given by 

02,0 (y) = D 2 , y J + 

where we define for 0(e m ) in the expansion 



n=l 



y 2 (Ua2 + U2h - u 2 a - 3D 2 , /i 2 





2h 








Jo dy 




f 27T d^ m 


2^ 


Jo dy 



yu2o, 



=o dx, 



\y = h 



(31) 

(32) 

(33a) 
(33b) 



as the mean components of the horizontal boundary conditions. The vertical boundary conditions cannot have a mean 
component and therefore do not contribute to the zeroth Fourier mode of the solution. 

Again to resolve the unknown constant we turn to dynamical considerations. Utilizing Eq. (j2"8"]l and evaluating at 
order e 2 we obtain 



«ao(0) = 4oW, 
therefore Z?2,o = 0. The force on the top sheet is then given as 

f2x = (""20 _ U 2/i - UA2)- 

h 

The mean components of the horizontal boundary conditions must then be evaluated, the lower 



U2Q = h 



2jt 



dy 2 



lv=o +2 



1 



sf(*Y ^ I g'itydt 



2u 



dx. 



The term in the square brackets clearly integrates to zero hence we are left with 



u 2 o 



1 

2Wo 



9{x) 



dy 2 



=o dx, 



which, using orthogonality of Fourier modes, gives 

n=l 

Similarly U2h is given by 



U 2 h 



1 

'2^ 



g(x + 



dy 2 



\y = h dx, 



(34) 



(35) 



(36) 



(37) 



(38) 



(39) 



which may be evaluated to give 



U2h 



1 ^ 

-J2\ («»<„ (h) + Pnb'{^ n {h)) cos(r 



+(/3 n a" n (h) - a n b'i {h))sm{nq 



(40) 
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Further, by considering P n and Q n , given by Eq. (|22[) and Eq. (l23l) respectively, and observing that 2P^(h) = — QJ'(O) 
and 2P£(0) = Q"(/i) it can be shown that 



(0)+«„ (0) 



1 a" „ (h) + f3 n b" „ (ft) cos(?k; 



P n a" n (h) - a n b" n (h) sm(n<, 



(a 2 n + Pl) 2P%(h) ~ Q'^O) - 008(^) 2^(0) - Q^(ft) 



0. 



(41) 



for all n. 

The force on the top sheet is then equal to 



27T rr 

J2x = -^rUA2- 

h 



(42) 



Here again we see that when we allow the swimmers to move in a force free manner then U A2 = and hence there is 
no synchronization at 0(e 2 ). Note that we have not specified the Fourier coefficients of the of the waveform, and this 
result is therefore valid for any waveform g(x). 

Sadly, due to the e — > — e symmetry of the model there cannot be any force at 0(e 3 ), and therefore we expect the 
force to arise at best at 0(e 4 ). 



3. Flow at 0(e 3 ) 



The third order component of Eq. @ is 



With the third order boundary conditions 



3 \y=0 



vVs = o. 



dy 



2 ' V dy 2 

Jl„\2 r 



l»=0 



2tt 







W> 3 | y= s = e y C/ A 3 - g(x + 
g'{x + <P 



8^ 2 \ gjx + t) (d 2 4> 
dy J + 2 {dy 2 



2tt 



9'(x + cf>) 2 - — I g'(x + cf>) 2 dx 



The solution can shown to be of the form 

oo 

ip3 = a 3 ,o{y) + X! [ a 3,"(y) cos(nx) + b 3 , n (y) sin(na;) 

where 



03,0(2/) = D 3i o2/ 3 + 



(Ua3 + U3h ~ u 30 - 3D 3 , h 2 ) y 2 
2h 



\y=h 



yu 3 o- 



(43) 



(44a) 



(44b) 



(45) 



(46) 



Again we must resort to dynamical considerations to solve for D30 and Ua3- From Eq. ([28]) we find a' 3 ' (0) = a' 3 ' (h) 
and hence D3.0 = 0. The force again takes the form 



2tt 

hx = ~y~ (W30 - u 3h -U A 3) 
h 



(47) 



If the swimmers are force free we see U a?, = U30 — but due to the e — > — e symmetry of the geometry we must 
have Ua3 — (in the example we consider in Sec. IVIll u^n = v,3h = 0)- 
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4. Flow at 0(e 4 ) 



The fourth order component of Eq. © is 



V 4 V>4 = 0. 



(48) 



The boundary conditions at fourth order are given by 



W>4 \y=0 



-g( x ) 



1 dy J 2 \ dy 2 J 6 V cty 3 



|j/=0 



V-04 |y=s = -5(a; + 0) 

<?'0 + 



9^3 



V 



2 ff ^ 



2 \ dy 2 J 6 



2tt 



27 / 



-e !/ ?7A4 



16tt 7 



2jt 



The solution is of the form 



tpi = 04,0(2/) + X! [ fl 4,"(y) cos(nx) + bi, n (y) sin(nx) 



(49a) 



(49b) 



(50) 



However we are only interested in its mean component, i.e. 



040 (y) = D4,oZ/ 3 + 



V 2 (Uaa - 3D 4 , /i 2 + u 4h - u 40 ) 



2h 



y«4o, 



From Eq. (|28p we obtain 040(0) = a'l (h) and hence D4.0 = 0. The force on the upper sheet is then 

2ir 

Ux = -T~(U40 - U4h - U A4 ). 

h 



(51) 



(52) 



Setting Ua = gives rise to a phase-locking force in the static case, /J (we use the superscript s to avoid confusion), 
given by 



2-7T 

h 

For free-swimming we thus see that the relative swimming sped is given by 

Noting that d<f>/dt = we therefore get an equation for the time-evolution of the phase as 



(53) 



(54) 



J t =-^f s x =e\u ih -u i0 ) + O(e% 



(55) 



Importantly, the formulae for U40 and Uih, defined in Eq. (|33[) . are too unwieldy for the most enterprising appendix even 
for simple g(x), and hence are not stated explicitly (although straightforward to obtain with a symbolic calculation 
package). In Sec. IVIII these analytical results for both the phase locking force and the dynamical problem will be 
compared with a full numerical solution using the boundary integral formulation. 
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C. Energy dissipation 

The energy dissipation rate between two sinusoidal sheets was originally computed by Taylor at leading order in 
the wave amplitude [ljj]. Here we restate his results for a general traveling wave. The energy dissipation per unit 
width in the fluid is equal to the rate of work of the sheets against the fluid 



E 



/ (u • er • n) \ y=yi dS - / (u ■ a • n) \ y=y . 2 dS. (56) 



Expanding the integral in e we find to leading order 

E = e 2 J J ' g'( x )(- pi + 2 ^j \ y=0 dx 



fj g'(x + 4>) (- Pl + \ y=h dx. (57) 



Expressing the pressure in terms of the stream function and integrating by parts yields 

E=-e 2 J g(x)-^- \ y=0 dx + e 2 J g(x + <j))-^-± \ y= - h dx. (58) 

We already know the form of these integrals (indeed they are equal) from the analysis of the force at 0(e 2 ), and we 
find 

oo 

e = ^y^^i+pd 

n=l 

x [Q«'(0) - 2P':>(h) cos(n0) (Q'i'(h) - 2P^"(0)) ] , (59) 

which we can evaluate to get 



E = 2ne 2 n 3 (al + (3 2 n ) [A{nh) - cos(n(f)B(nh)] , (60) 

where 



n=l 



MO = (61a) 



It, + sinh 2^ 
sinh 2 £, - i 2 
2^ cosh£ + 2 sinh ^ 
sinh 2 ^ — £ 2 

Setting Pi = 1 and all other coefficients to zero in the above yields Taylor's result for pure sinewaves [l^ |. 

In the limit £ — > oo, we see that A — > 2, B — > 0, and the ratio B/A decays exponentially. This tells us what we 
intuitively expect: When h is large, the phase difference has little effect on the energy dissipation, and also the phase 
difference has a weaker effect on the energy dissipated by higher Fourier modes. Conversely, we expect that when the 
separation is small, the phase angle would have an large influence on the rate of working of the swimmers, and indeed 
when £ — > 0, we zee B/A — >• 1 as both A,B~^ 12£~ 3 (keeping in mind that we have implicitly assumed e ^ h). 

Importantly, because A and B are both positive and monotonically decaying functions with £, we know that in- 
phasc swimming, cf> = 0, is a global minimum for the energy dissipated in the fluid. In addition, given that we have 
the symmetry g(x + ir) = —g(x), this restricts us to odd Fourier modes, and thus the out-of-phase configuration, 
4> = it, is a global maximum. Taylor's dissipation argument fl4| extends thus to arbitrary waveforms. 
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V. LUBRICATION LIMIT 

A second insightful limit to consider is the one in which the sheets are so close together that their mean separation 
is much smaller than the wavelength of the oscillations, h <C A. In this lubrication limit the Stokes equations are 
substantially simplified, permitting analytical solutions. The main results of this section were previously summarized 



in a letter by the authors 26 1. 



A. Lubrication equations 

In order to facilitate this limit we must rescale the governing equations. We nondimensionalize vertical distances 
by y* = y/h, and horizontal distances z* = kz, while assuming that 5 = kh <C 1. The instantaneous position of the 
sheets is therefore given by y* = a*g(x*) and y% = l + a*g(x* + </>), where a* = a/h and again x* = z* — t* is the wave 
variable. Nondimensionalizing the continuity equation we find that if the horizontal velocity is given by u = cu* then 
the vertical velocity must be v = Scv*. The Stokes equations then yield the lubrication equations to leading order in 
<5: 

du* _ dv* . 
dx* dy* ' 



djf 
dy* 
dp* 



0, (63) 
ft 2 ,,* 

(64) 



dx* dy* 2 ' 

where p* = 8 2 p/ /MjJ. Forces (per unit depth) are nondimensionalized as /* = fS/fic, while energy dissipation rate per 
unit depth is E* = S 2 E/ [lujch. We note that if g approaches a singular geometry we would leave the realm of validity 



of the lubrication approximation [34l. I35| . We now drop the * notation for convenience. 

We look to solve this problem in a frame moving with the wave speed of the bottom sheet. The boundary conditions 
in the lubrication limit are then given by 

u(x,y = ?/i ) = -1, (65a) 

v(x,y = yi) = -y[, (65b) 

u(x, y = y 2 ) = Ua - 1, (65c) 

v{x,y = j/2 ) = -y' 2 . (65d) 

Hence we see that in the lubrication limit the boundary conditions are identical to those of an extensible sheet. 
The full problem, regardless of whether extensible or inextensible boundary conditions are used, will collapse to the 
following lubrication results in the limit S <C 1. 

Given the above boundary conditions, the solution for the velocity field is found to be 

u{x,y) = \^-{y - y\){y - 2/2) + U A - — — - 1. (66) 
2dx 2/2 - yi 



If one integrates the continuity equation one finds 



in 



) dy = y! l -y , 1 , (67) 



which then gives 



dx dx 

If no relative motion, Ua — 0, then the flow rate Q between the sheets is constant 



d Q tt d y 2 faa\ 
Ua-t-- (68) 
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B. Hydrodynamic force 



We first characterize the force generated when Ua = in order to determine the location and nature of the fixed 
points for the phase difference between the swimmers. With Ua = then Q = const., and we find 



V2 a 1 d Pu3 , 

udy = n~]~< 1 ~ 



12 dx 



where h = 1J2 — yi ■ We now exploit the periodicity of the system to obtain the value of Q by noting that 



dp 
dx 



dx = -12 



1 

h 1 



We thus have 



where 



dx - 12Q 



h 

V 



1 

If' 



dx = 0. 



The pressure gradient is therefore given by 



h = 



£ = 12 

dx 



h •'dx. 



1 



h 3 I 3 h 2 

The force per unit depth on the upper sheet is given by 



fx = e x - / a ■ n \ y=y2 dS, 
Js 



where the curve S is defined by y = 1/2 ■ Evaluating Eq. (|74|) gives 

r 2ir 



fx 



(I 



dx \ ox J \oy ox 



\y=V2 dx. 



We keep only the O(l) terms in the lubrication limit 8 <C 1 which yields 

Jo v dx d v 

Exploiting the periodicity of the problem through integration by parts 

dp du 
' dx dy 



361 allows us to recast the force as 



»r-5: lw=to dx - 



(69) 



(70) 



(71) 



(72) 



(73) 



(74) 



(75) 



(76) 



(77) 



Substituting in Eq. (|66j) and Eq. ()73|) . and noting any constant multiplying the pressure gradient may be discarded, 
we find the force to be given by 



f x = 6a J (j^ - J^j [d( x + 4>) + 9{x)] dx. 



(78) 



C. Fixed points 

By symmetry, we found earlier that there are always fixed points at (f) = 0, tt. This is easily confirmed by evaluating 
Eq. (f78|) . For cj) = 0, h is constant, and thus ^/h^Iz — l/ti 2 = 0, leading to f x = 0; for <f> = ir, we have g(x+ir)+g(x) = 
by symmetry, and again f x = 0. 
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In order to determine their stability, we can expand the force, Eq. (l78l) . about these fixed points. Letting (j> = No- 
where <p' <C 1 we obtain near the in-phase fixed point, 4>q = 0, 



f X0 = -72aV 3 r g{x)g'{xfdx + 0{4> li ). 



II 



In contrast, near the opposite-phase fixed point, 4>q = n, we get 

[*" 9'{x? ( 

/o {l-2ag{x)Y \(l-2ag(x)) J 3 

where we have defined 



f x „ = 6aV 3 9 ' {x) ', d ( - * , J -± ) dx + O(0' 4 ), 



[■2-ir 

J n = (1 - 2ag(x))- n dx. (81) 
Jo 

If we then assume a <C 1 then Eq. reduces to 

/^«72aV 3 / 5 W(^) 3 ^ + O(0' 4 ). (82) 
Jo 

We see then that for small amplitude waves, and small deviations in phase about the fixed points, the force about 
the in-phase configuration ((/> = 0) is equal and opposite to the force about the out-of-phase configuration (<j) = n). 
Unless both of them are neutrally stable (which is the case if the waveforms are too symmetric, see Sec. we 
therefore obtain the important result that, for a given waveform, one fixed point will always be stable, while the other 
one will always be unstable. To determine which one is the stable point, one has to evaluate the geometric integral, 
A = J Q 2lr gg' 3 dx. If A < then the fixed point at <fi = is stable, while it is the one at = 7r in the case A > 0. Stable 
passive hydrodynamic synchronization thus always takes place for swimmers with asymmetric waveforms. 

As a side note, we can also expand the force (|78|) in powers of a <C 1, and we see that the leading order contribution 
is fourth order in amplitude, given for general </> as 

U a -36a 4 f (g{x + <£) + g(x)) (g(x + cf>) - g{x)) * dx, (83) 



plus terms at 0(a 6 ). We see that in the small amplitude limit there are only two fixed points for nontrivial waveforms 
g. The fourth-order scaling of the hydrodynamic force, Eqs. ([79| . ([82]) and (|83|) . is reminiscent of the small-amplitude 
calculations from Sec. II VI showing that no force can occur at second order in the wave amplitude, but a nonzero force 
does come at fourth order. 



D. Energy Dissipation 

The energy dissipated by viscous stress in the volume V of fluid between the sheets by is given by 

E = I cr : Vu dV. (84) 
Jv 



In the lubrication limit, assuming unit width, the energy dissipation over one wavelength is then 

/" 27r f V2 fdu 



E = r f 2 (ir) d V dx > ( 85 ) 



and given Eq. (|66p we have 



E=^r^( d ±\ dx, 



12 ./o V dx 
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which is explicitly 



(87) 



We see the energy dissipation is non-negative and identically zero when = (i.e. when h is constant) and hence 
this must be a global minimum. 

If we let = 0o + 0' where 0' <C 1, we find near the in- phase conformation, O = 

E Q = l2a 2 (t)' 2 g'(x) 2 dx + 0((j)' 4 ), (88) 
Jo 

and the energy increases quadratically with the slope of the wave from zero when = 0. Near the opposite-phase 
conformation, 0o — 7T, we get 



E v = 12 Ji 



72 
J 3 



120 



g'{xf 



1 



G 



o | (1 - 2ag{x)f 
h ( 1 J 2 



(1 - 2ag(x)) J 3 V(l -2ap(a;)) J 3 
If we further assume that a <C 1 we see that 



dx } + O(0' 4 ). 



E„ « 12a' 



2tt 



[4 ff (x) 2 - 3 '(x) 2 0' 2 ]dz, 



(89) 



(90) 



hence for any waveform g{x) the energy dissipated between the sheets is maximum in the opposite-phase conformation, 

= 71". 

Finally, if we expand the energy dissipation, Eq. (|87j). in small amplitude for general 0, we find 

[g{x + <j>) ~g{x)] 2 dx + 0(a 3 ). 



12a' 



(91) 



We can see clearly again that the energy dissipation is a global minimum when = and maximum when = tt due 
to the g(x + tt) — > —g(x) symmetry of the waveform; this is in agreement with the previous small amplitude results 
for arbitrary separation. 

An important consequence of the previous results is that, although the nature of the fixed points depends on the 
swimmer waveform, the location of the minimum of energy dissipation does not. The conformation of minimum 
energy dissipation is not necessarily stable: depending on the waveform geometry, the opposite-phase conformation, 
= 7T, may be stable (specifically, when A > 0) yet it is the one corresponding to a maximum of energy dissipation. 

Experimental evidence shows that spermatozoa cells synchronize at the in-phase conformation (and indeed A < 
for the linearly increasing sine waves indicated by Rikmcnspoel [31(). However, we find at least one instance, in the 
figures in Ref. [12| , which show spermatozoa flagella seemingly synchronized in opposite-phase (although no mention 
of phase difference is reported in the text). 



E. Dynamics 

After calculating the hydrodynamic force, we now look to solve for the time-evolution of the phase. We thus assume 
that the sheets are force free, f x = 0, and find the corresponding value of From Eq. (|68j) we know 



*r*- aJ s- (92) 
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Integrating in x and evaluating the integral in y gives an expression for the pressure gradient as 

dp 6U A - 12 l2U A y 2 + C 



dx 



h 2 



h 3 



(93) 



where C is a constant of integration. We find this constant by exploiting the periodicity of the pressure field, leading 
to 



C = (&U A {h-2K)-l2I^/h 



(94) 



where K = J Q w y 2 h 3 dx. 

The force on the upper sheet is given by 

fx 



I (lt iy2+m) -^) dx - 



(95) 



We then solve for U A hy enforcing that the sheets are force free. It is worth noting that when we set U A = 0, 
we retrieve the force from the static case given by Eq. (|78p . which we label here /J to avoid confusion. Now since 
U A = —d<f>/dt we find that the phase evolves in time according to 



where the mobility, A4 , is given by 



M- 1 = 



dt 



1 1 



2y 2 + 



In - 2K 



h 



3(j/2 + yi) } x. 



(96) 



(97) 



As physically expected, the rate at which the phase changes, Eq. (|M|) . is proportional to the static force, /J, which 
would be applied if the sheets where not permitted to move. The result is a first-order integro-differential equation 
for 4>. 

Expanding Eq. (|96[) for small amplitude, n«l, we find 



d(j> 36a 4 
~dt 2tT 



(g(x + <t>) + g(xj\ (g(x + 4>) - g{x) 



dx, 



(98) 



plus terms at order a 6 . We thus see that d<p/dt ~ — j s x j2-K for small amplitude, and hence the mobility becomes 1/2-k 
in this limit. Notably, the result given by Eq. ([M]) is the same as the one given by Eq. ([53]) after proper dimensional 
rescaling. 

b' and obtain 



We now expand near the fixed points by letting tfi = (fio 

36 



dt 



7T 



(99) 



with a positive sign for 4>o — 0, and negative cf>Q = it. Solving this differential equation gives the exact phase dynamics 
for small times as 



sgn($) 



=F 72a 4 At/Tr 



(100) 



where <p'(t = 0) = fy. In the case of a stable fixed point, we thus get that the typical time for synchronization scales 
as t ~ 1/a 4 1 A| , and thus the phase- locked state is reached faster for waves of larger amplitude (a increases), and 
larger asymmetry (\A\ increases). Note that, as a difference, the typical time for synchronization in a viscoclastic 
fluids scales as the inverse square of the wave amplitude 30 1 . 
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VI. BOUNDARY INTEGRAL FORMULATION 

The boundary integral method may be used to address numerically the synchronization between the swimmers for 
shapes of arbitrary amplitude, as well as confirm our asymptotic results. We present in this section the principle of 
the method and our implementation of it, which is quite similar to that given by Pozrikidis in his study of peristaltic 
pumping [371 ] . and hence will be brief. The equations in the section are nondimensionalized similarly to what was 
done in Sec. UH 

Consider any two solutions to the Stokes equations, {u, cr} and j_u, er} with no associated body forces for any closed 
surface S of outward normal n. The Lorentz reciprocal theorem [381 ] gives the equality 

/ (u-a-u-(r)-ndS = 0. (101) 
Js 

If we take for u and cr in Eq. (|101l) the fundamental solutions for two-dimensional Stokes flow 

u(x) = i-G(x) • f (xo), (102) 
<r(x) = i-T(x) • f (x ), (103) 

for the velocity and the stress at the field point x, due to the point force f at xo, where x = x — xo and the 
two-dimensional Stokeslet G and stresslet T are given by 

V V 

G = -Iln(|x|) + — (104) 
l x l 

T = -4^f , (105) 
l x l 

then one obtains the boundary integral formulation of the two-dimensional Stokes equations for the velocity field 
within the fluid domain, V , and on the boundary, S, respectively, 

u(x )| xoe y = — — / (u(x) • T(x) ■ n - b(x) • G(x)) dS(x), (106) 
J s 

u(x )| xoes = y f (u(x) • T(x) • n - b(x) • G(x)) dS(x), (107) 
J s 

where we have used b = cr ■ n. 

Since the problem we consider is 27r-pcriodic, we can reduce the domain of integration S to a single period by using 
an infinite sum of periodically placed Stokeslets and stresslcts, 

G* = -iMWJ + S. (108) 

n=-oo l X ™l 
oo ~ 

TP = J2 -4 X " X " X " , (109) 



where x„ = {xq + 27rn, yo} f^ . As shown in Ref. [3] these may be expressed in closed form by using the summation 
formula 

B= J2 ln(|x„|) = -ln[2cosh(£ )-2cos(£ )]. (110) 
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We can then construct the elements of G p and T p using B and its derivatives as follows: 



G p 

^ XX 
QP 

W X y 

QP 

yy 

rpp 
rpp 

xxy 

rpp 

xyy 

rpp 

yyy 



-B - d v B 



1, 



yd x B, 

-B + yd y B, 
-2(2d x B + yd xy B), 



-2(d y B - 
2yd xy B, 
-2(d y B - 



ydyyB), 



ydyyB), 



(111) 



where the Stokeslet and stresslet are invariant under permutation of its indices 

Following the approach outlined by Higdon [4l|, the boundary S (the surface of each sheet as the sides cancel) is 
discretized into 2N straight line elements S n . We assume the stress b and the velocity u are linear functions over a 
particular interval (b — > b„, u — > u„) and then collocate xo at each of the 2N segments, Xo — > x m , to obtain a system 
of N equations with TV unknowns. The periodic Stokeslet and stresslet are regularized by subtracting off from them 
their non-periodic counterparts and then adding back the difference; the two-dimensional Stokeslet and stresslet are 
then integrated analytically. We hence have 



u(x m ) = 



2N 

— y 

71 = 1 



u„ • (TP - T) • n n dS n 



u„ ■ T • n n dS„ 



h n ■ (GP - G) dS n 



GdS n 



(112) 



The regularized integrals have a removable singularity at x = x m which is obtained by Taylor expansion. The 
boundary integral formulation is thereby reduced to a linear system that can be inverted using standard techniques 
to obtain the stress b. The force on the top sheet is then given by integrating the stress 



2iV 

E 

! = JV + 1 



\) n dS n 



(113) 



In order to solve for the dynamical problem we let the boundary condition be represented by a sum of a surface 
deformations and an unknown rigid body motion, u->u n + Ua^x, on the upper sheet. The additional unknown, 
is found by enforcing that the sheets are force free, f x = 0. 

The numerical procedure was validated through convergence tests and by reproducing previous results for shear 
flow over sinusoidal surfaces 40 1 . 



VII. RESULTS 

We now present in this section the results of both our asymptotic and numerical calculations to address the 
synchronization of specific waveforms. For illustrative purposes we restrict ourselves here to the family of waveforms 
described by the function, 

g(x) = sin.T + 7cos3o;, (114) 

i.e. Pi = 1, az = 7, and all other modes equal to zero, as illustrated in Fig. [4] In essence these shapes are geometric 
perturbations (small for small 7) to Taylor's sinusoidal swimming sheet. They have a broken front-back symmetry 
when 7 is nonzero. Reversing the sign of 7 is equivalent to reflecting the geometry of each wave about the vertical 
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axis, {7 — > —7} = {x — > —x + 7r}, which itself is equivalent to reversing the kinematics of the problem. In other 
words, changing the sign of 7 changes the sign of the forces on the sheets which leads to stable fixed points becoming 
unstable, and vice versa. In addition, the simple form of g(x) allows us to obtain some explicit formulae from the 
general theory presented in Sees. HVIandlVl 

In the lubrication limit, the geometric parameter A = gg'^dx = —2tt^/ controls the evolution of the phase near 
fixed points. We see that 7 > gives A < 0, which leads a stable fixed point at = and unstable at = it. By 
symmetry, 7 < necessarily gives A > 0, and thereby exchanges the location of the stable and unstable fixed points. 
In addition, from Eq. (|83[) , we have that when a < 1 the phase locking force is given by 

f x = 1447ra 4 7 sin 3 0, (115) 

which is linear in the asymmetry and quartic in the amplitude, and leads to a time-evolution of the phase as 

— = -72a 4 7 sin 3 0. (116) 
at 

The energy dissipation in the lubrication limit, for a <C 1, Eq. (pTT|) . is 

E « 247ra 2 [l - cos0 + 7 2 (l - cos 30)] . (117) 
Similarly, in the small amplitude limit, Eq. (|60[) yields 



E « 27re z 



A(h) - B(h) cos + 3 3 j 2 (A(3h) - B(3h) cos 30) , (118) 



where the functions A and B, given by Eq. (161[) . introduce a dependance on the separation h, and Eq. (|118[) reduces 
to Eq. (|117j) when h is small (after accounting for the separate scalings). 

We see clearly that the energy dissipation rate is invariant under 7 — > —7 and further, when we are assuming that 
7 is a small, the change in the energy dissipation due to the asymmetry is also small, 0("y 2 ). 



A. Comparison between asymptotic and numerical methods 

In the small amplitude limit pV[) we have explicitly assumed that e C 1, and also implicitly that e <C h. The 
lubrication limit (|V| effectively captures the physics of the problem when the sheets are quite close together, i.e. the 
limit h 2t:. If we want in addition the phase, 0, to be able to span the range of all possible values then we also get 
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FIG. 5: (Color online). Force, f x , vs. phase difference, <j), for an asymmetry of 7 = 0.1, in the lubrication limit (dashed 
lines, top figure only), small amplitude limit (solid lines, both figures), and using the boundary integral method (symbols, both 
figures). Top: fixed amplitude, e = 0.01, and varying swimmer-swimmer distances, h; blue circles: f x , h — 0.2; red squares: 
10/3,, h = 0.4; green diamonds: fx/e, h = 0.6; gray triangles: 2f x /e, h = 0.8. Bottom: Fixed separation distance h = 1 and 
varying waveform amplitudes; green circles: Tvf x /e, e — 0.1; red squares: nf x with e = 0.2; blue diamonds: f x , e = 0.4. We 
observe a gradual breakdown of the lubrication approximation for increased separation, h (top), and of the small amplitude 
expansion for increased amplitude e (bottom). Note that forces have been scaled for display purposes. 



the geometrical constraint e < h/2 (or, in terms of lubrication variables, a < 1/2). There exists therefore a regime in 
which both asymptotic approaches are valid, namely the limit £ < /i « 1. 

As a validation of our methods we plot the analytical results from both asymptotic limits, together with the 
numerical results, for such a regime in Fig. [5] (top). The static force on the top sheet, /J, is shown for the waveform 
of Eq. (|114j) with an asymmetry of 7 = 0.1 and wave amplitude e = 0.01 (in this plot the forces have been scaled for 
display purposes only, see caption). The solid lines represent the small amplitude limit, dashed lines the lubrication 
limit, and symbols are for the numerical data obtained from the boundary integral method. The results from all three 
methods agree quantitatively for small swimmer-swimmer separation, h. As the value of h increases, the lubrication 
results start to deviate, but the small amplitude results remains accurate (recall that e <C 1 in all cases). 

For larger values of the separation distance between the swimmers, the lubrication results cannot be applied, but 
the small-amplitude asymptotics, Eq. ([55]) , remain valid as long as the wave amplitude is small. The value of the 
static force is compared to the numerical results in Fig. [5] (bottom) for large separation, h = 1, and as function of the 
wave amplitude, e. The agreement between the two is excellent for e = 0.1, but they deviate quantitatively for larger 
wave amplitudes (although the correct order of magnitude, and dependence on </>, is obtained). 
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FIG. 6: (Color online). Dependence of the force, f x , on the phase cj> f° r varying separation ft with an asymmetry 7 = 0.1. 
Top: small dimensionless amplitude, e = 0.1. The solid lines are obtained in the small amplitude limit while symbols are for 
boundary integral computations; blue circles: ef x , ft = 1; red squares: 2f x , ft = 2; green diamonds: nf x /2e, ft — 3; purple 
down triangles: fx/e 2 , ft = 4; orange up triangles: fx/f 3 , ft = 5; gray stars: fx/ire 4 , ft = 6. Increasing the distance between 
the sheets introduces an additional fixed point not present in the lubrication limit, and its position moves with ft. Bottom: 
numerical results using the boundary integral method (solid line and symbols) in the case of high amplitude waves, e = 1. 
Black triangles: f x , ft = 3; red diamonds: 27r/ !E , ft = 4; green squares: lOrr/a;, ft = 5; blue circles: 100 /r, ft = 6. Forces have 
been scaled for display purposes. 

B. Stability 

When we introduce a variable separation between the swimmers, h, and thus go beyond the small h limit from the 
lubrication approximation, we get that the number of fixed points and their nature does not depend solely on the 
waveform geometry, but actually also on the swimmer-swimmer distance. In Fig. [6] (top) we show the dependence of 
the static force on the phase, for an amplitude e = 0.1 and an asymmetry 7 = 0.1, as we vary the separation between 
the swimmers h (line: small-amplitude asymptotics; symbols: boundary integral computations). A fixed point is a 
conformation with phase difference <j> such that f x (4>) = 0; if the slope of the force is positive the fixed point is stable, 
while a negative slope indicates an unstable fixed point. What we see in Fig. [B] (top) is that increasing the separation 
between the sheets from the small h values in the lubrication limit gives rise to an additional fixed point. In the case 
illustrated in Fig. [5] (top), this new fixed point is always unstable. It first appears near <p = ir (leading to the fixed 
point at 4> = 7T becoming stable), moves toward <f> = when the separation distance between the swimmers increases, 
and eventually merges with = 0, which then turns to an unstable point, at a critical value of h. 

In Fig. [7] we display the location of the fixed points explicitly as a function of h in the small amplitude limit for 
7 > 0. In this limit, the force is linear in the asymmetry, 7, therefore the location of the fixed points is invariant 
under a linear scaling of the asymmetry, 7 — > 67 where b > 0, while the nature of the fixed points changes with a 
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FIG. 7: (Color online). Location of the fixed points (value of the phase (j> such that fx{(j>) = 0) as a function of h as obtained 
in the small amplitude limit for 7 > 0. The solid blue lines indicate a stable fixed point whereas the dashed red line indicates 
an unstable fixed point. The position of the unstable fixed point moves with h: it is created for small values of h near <j> — n, 
migrates to the left, and merges with (f> = at the critical value h « 5.65. In the opposite case where 7 < 0, the stable fixed 
points become unstable and vice versa. 

change of the sign of 7. The appearance of a new fixed point, described in the previous paragraph, is apparent. As h 
tends asymptotically to zero, <fi — is stable (blue solid line) while = tt is unstable (red dashed line), which is the 
lubrication result. For intermediate values of h, both and tt are stable, and the new fixed point moves from tt to as 
h increases. It merges with <fi = for a critical distance between the swimmers (h « 5.65 for our choice of waveform), 
at which point = becomes unstable, and remains so for larger values of h. As expected, upon a reversing the sign 
of 7, stable fixed points become unstable and vice versa. 

Further analysis of the equations of motion shows that the additional fixed point that arises when h is past the 
lubrication limit is a direct consequence of the inextensible boundary conditions. In the lubrication limit, the boundary 
conditions are extensible insofar as the there is only a vertical component, however away from this limit there arises 
horizontal motion to maintain inextensibility, and it is precisely this horizontal motion which leads to the additional 
dynamic complexity. Conversely for extensible boundary conditions, Eq. (|10[) . the fixed points remain unchanged 
from those in the lubrication limit. 

Using the boundary integral formulation it is possible to extend these results to large amplitude waves. In Fig. [6] 
(bottom) we show the horizontal force on the upper sheet as a function of the phase between the swimmers for various 
mean swimmer- swimmer separation but now with e = 1. The results are qualitatively similar to those obtained in the 
small amplitude limit, with the occurrence of a new fixed point, unstable, and moving from <f> = tt to </> = as the 
separation increases. A difference we do observe between small and large amplitude is that the location of the fixed 
point is no longer invariant under a change in the asymmetry factor, 7. In Fig. [8] (top) we show that an increase in 
the asymmetry factor leads to a small, but nonzero, migration of the mobile fixed point toward tt. A similar drift is 
obtained with an increase in the waveform amplitude (Fig. [8j bottom) . 



C. Dynamics 

The time-evolution of the phase is given in general by the integro-differcntial equation 

^ = -M(4>)f^). (119) 

As noted above, in the small amplitude limit the mobility becomes independent of the phase, M. = h/2ir, hence in 
that case the dynamics is completely set by the static force. Note that the mobility is never zero so no additional 
fixed points arise from Eq. (j!19[) . 



24 




FIG. 8: (Color online). Plot of the force on the top swimmer, /„ as a function of the phase difference, <f>, using the boundary 
integral method with h — 4. Top: e = 1 for varying asymmetry; blue circles: 7 = 0.05; green squares: 7 = 0.1; red diamonds: 
7 = 0.2; black triangles 7 = 0.3. We see that for large amplitude waves the force is no longer linear with asymmetry as 
evidenced by the moving of the middle fixed point. Bottom: 7 = 0.1 for varying large amplitude waves; blue diamonds: 107r/ x , 
e — 0.5; green squares: 2nf x , e = 1; red circles: f x , e = 1.5. We see that for e < 1 the location of the middle fixed point remains 
close to the small amplitude limit, while it has drifted significantly for e = 1.5. Forces have been scaled for display purposes. 



In the lubrication limit we know that there exist only two fixed points, and the location of the stable fixed point 
depends only on the waveform asymmetry. In Fig. [5] we plot the time-evolution of the phase in this limit. We see 
that if the system is symmetric (7 = 0), indicated by the black solid line, then the phase remains constant in time. 
This corresponds to the no-synchronization situation discussed in Sec. IIIIl When we introduce an asymmetry, 7 7^ 0, 




FIG. 9: (Color online). Time-evolution of the phase, from <f>o = 7r/2, in the lubrication limit with 8 = 0.1. The dashed line 
indicates 7 = while the solid line indicates 7 = ±0.1, a = 1/4. The dashed line has 7 = 0.05, a = 1/4 and the dash-dot line 
has 7 = 0.1, a = 1/8. 
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FIG. 10: (Color online). Time-evolution of the phase, from the initial condition 4>o — 37r/4, for h = [1 — ► 4] in the small 
amplitude limit with e = 0.1. A solid line indicates 7 = 0.1 while a dashed line indicates 7 = 0.2. 

then the two swimmers phase lock over time. When 7 > then A < and the system evolves to a stable in-phase 
conformation, and opposite-phase for the converse. Given that the amplitude, a, is reasonably small for all curves 
(we have the geometrical constraint a < 1/2), we observe roughly the same dcpcndancc of the typical time scale for 
phase locking, t, on the wave asymmetry and amplitude as in the small-amplitude limit (for which t ~ o — 4/ y — 1 , see 
Eq.[99]). 

We have seen above that with an increase in h comes additional fixed point, and thus we expect the phase dynamics 
to depend similarly on h. In Fig. 1101 we plot the time-evolution of the phase in the small amplitude limit for various 
values of the swimmer-swimmer distance in the case where 7 > 0. Given that the phase mobility, M, is independent 
of the asymmetry, and that the force is linear in 7, we find that the time scale for synchronization scales with the 
inverse of the asymmetry factor, i.e. t ~ 7 , as it does when h -C 1. The final stable swimmer-swimmer conformation 
can be understood simply by recalling the force plot in Fig. [JJ If the initial conformation is to the left of the moving 
unstable point, then the sheets evolve to </> = 0, while they start to the right they evolve to <f> = it. If we reverse 
the asymmetry of the waveforms, 7 < 0, then the converse is true, the fixed point which varies with separation 
represents the only stable conformation for intermediate values of h and we obtain synchronization to a fixed finite 
phase difference, < (j> < it, as is observed in the metachronal beating of cilia. 

A similar plot is shown in Fig. llll in the case of large amplitude waves, using the boundary integral method, starting 
from an initial relative phase of <p = 37r/4 and with a positive asymmetry, 7 > 0. Again the essential physics is well 
captured by the small amplitude expansion: there exists a critical swimmer-swimmer separation below which the 
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FIG. 11: (Color online). Time-evolution of the phase, from 4>o = 37r/4, for large amplitude waves using the boundary integral 
method. The blue circles indicate h = 4, 7 = 0.1, e = 1.5; green squares: h = 4, 7 = 0.1, e = 1; red diamonds: h = 4, 7 = 0.05, 
e = 1 and black triangles: h — 2.4, 7 = 0.1 and e — 1. 
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FIG. 12: (Color online). Illustration of the flow field during synchronization for e = 1, 7 = 0.1 and h — 4 using boundary 
integral computations; (a) flow vorticity; (b) iso values of |u| 2 . Darker regions correspond to higher velocity and vorticity, and 
arrows indicate instantaneous velocity vector field. The two sheets start from the initial condition cf>o = 3n/4 (top panel). Time 
increases top to bottom with the phase <f) being equal to {3n/4, it/2, 7r/4, 0} in the four panels. 

sheets evolve to the in-phase conformation. This is seen in Fig. Qj] where with h = 2.4, e = 1 and 7 = 0.1 the sheets 
evolve to <fr = (in phase) whereas when the distance is increased to ft, = 4 the sheets evolve to <f> = tt (opposite 
phase). A waveform with a larger amplitude, e = 1.5, leads to a faster evolution of the phase than for e = 1 for equal 
asymmetry (7 = 0.1), which in turn evolves faster than for equal amplitude, e = 1, but lower asymmetry 7 = 0.05. 
We note however that for large amplitude waves, the effect of increasing the amplitude on the rate of phase change 
is drastically reduced; in the small-amplitude limit the rate of evolution is quartic with the wave amplitude and here 
we sec an effect which is less than cubic. Despite the reduction, the effect of amplitude is still strong and we observe 
drastically faster synchronization for order one amplitudes. 

To illustrate the flow driving the synchronizing dynamics shown in Fig. 111! we produce snapshots of the flow field 
for the h = 2.4, e = 1, 7 = 0.1 conformation in Fig. [12] We display the out-of-plane vorticity, to, in Fig. [T2k and the 
squared velocity field, |u| 2 , in Fig. IT2"b . Both plots are overlaid with arrows indicating the velocity vector field. The 
darker regions indicate higher vorticity and velocity in each plot respectively. Time increases from top to bottom, and 
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we show the instances where the phase angle is given by <j> = {37r/4, tt/2, tt/4, 0}, corresponding to relative velocity of 
the sheets, U A = {0.0325,0.1223,0.1097,0}. 



VIII. DISCUSSION 
A. Summary of results 

In this paper we have considered, as a model for the synchronization of flagellated cells, the relative motion of 
two free-swimming planar parallel sheets propagating waves of transverse displacement. We showed that due to 
the kinematic reversibility of the Stokes equations, waveform conformations with both vertical and horizontal axis 
symmetry could not yield synchronizing dynamics. When we break vertical axis symmetry, the phase of the system 
evolves to stable dynamic equilibria whose position is entirely dependent on the geometry of flagellar waveforms, and 
the distance between them. 

When the swimmers are close together we are able to make use of the lubrication equations and find two fixed 
points, in-phase and opposite-phase. The location of the stable point depends on the nature of the asymmetry, which 
is measured by an integral over the waveform geometry. If the front-back asymmetry of the geometry is reversed 
it is equivalent to reversing the kinematics of the problem which yields a reversal of forces. In other words, stable 
points become unstable and vice versa. In contrast, the energy dissipation is always a minimum for the in-phase 
conformation, and indicates the possibility of phase-locking into a conformation of maximum energy dissipation. 

An expansion in small amplitude is utilized to introduce inextensible boundary conditions and order one distances 
between the organisms. In this case there arise additional fixed points, whose location and nature depends on the cells 
geometry and separation. Among the possibilities is synchronization at a stable intermediate phase between in-phase 
and opposite-phase. 

Finally, we presented numerical results for large-amplitude waves using the boundary integral method. The compu- 
tational results indicate that between the lubrication limit and the small amplitude expansion, all the relevant physics 
can be captured analytically. However, since the phase locking force depends so strongly on waveform amplitude, we 
observe much faster synchronization for large amplitudes, as might be expected. 



B. Two-dimensional modeling and collective locomotion 

The two dimensional model used here is admittedly too simple to provide quantitative predictions for real microor- 
ganisms. However the simplicity allows analytic formulae to be derived and a mathematical structure elucidating 
the interaction between the bodies due to the Stokesian flow to be obtained, and gives an explicit description of the 
effect of symmetry breaking. The intuition garnered here may then be useful for more complex models, with finite 
three-dimensional bodies, that must be solved entirely numerically. 

We first note that, as a result of the two-dimensional approach, the viscous mobility of the cells in the direction 
perpendicular to that of the wave propagation is strictly equal to zero. For real microorganisms however this is not the 
case and hence fluid forces will determine the separation between the swimmers dynamically. Since swimming cells 
are force-free, the far-field velocity is typically a force dipole. In particular spermatozoa have a positive force dipolc 



(so-called pushers ll|). Far- field interactions between pushers tend to both align the major axis of pushers drive 
them together. Accordingly experimental evidence suggests that as spermatozoa synchronize they come together very 
tightly [5| (see also Fig. [IJ. In light of this, the lubrication limit is perhaps the relevant physical limit to consider for 
the phase locking of swimming cells. In contrast, cukaryotic cilia are attached to a substrate at a fixed distance which 
is varies depending on the organism flol |. 
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In addition, the two-dimensional limit offers one particularity, which is that the fluid between the swimmers (inner 
problem) does not communicate with that outside the swimmers (outer problem). The outer problem, that addressed 
by Taylor, is purely kinematic, in the sense that the speed of the sheet relative to the flow at infinity is uniquely 
determined without resorting to a balance of forces, unlike the inner problem. Further to this, because the outer 
problem cannot impose a force on the outer surface of the sheet, the forces are individually zero for the inner problem 
and therefore the inner problem (or even arbitrarily many inner problems) may be solved separately because a balance 
with the outer problem is not required. Now if a rigid body motion U = U&e x (due to the inner problem) is added 
to the surface deformations of the outer problem it has the sole effect of adding a plug flow solution to the swimming 
problem; in the zero Reynolds number limit a rigid body motion of two dimensional surface diffuses to infinity 
instantaneously. An interesting consequence of this is that when there arises a nonzero relative velocity, the idea of 
collective motion loses meaning, even for identical sheets, as in the frame moving with lower sheet (see Fig. [2]) we find 
different values for the flow at infinity, U when y — Y — oo and U + Ua when y — > oo. Further, if the sheets are different, 
then even if the inner problem demands J7a = 0, the outer problem on either side produces flows at infinity (in the 
frame moving with the sheets) which are distinct. However even in the three dimensional case, when the swimmers 
are the same, synchronization is clearly driven by the forces between the bodies and those forces will dominate when 
the cells are close; because of this we expect to garner the leading order behavior from analysis of the inner problem. 



C. Avenues for future work 

In the problem considered in this paper, we have prescribed a front-back asymmetry in the waveforms propagated 
in our model of flagellated cells in order to give rise to synchronization. Real eukaryotic flagella deform under applied 
(internal) forces, and this deformation may provide an additional mechanism of symmetry-breaking. Indeed recent 
experiments using rotating paddles suggest that elastic deformation is a key requirement to obtain synchrony for 



a geometry that is otherwise too symmetric to yield stable fixed points 42j. Flagella flexibility might thus be an 



important physical ingredient in the synchronizing dynamics, and we will address its relevance in future work. 
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